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Abstract 

The bird cherry-oat aphid (Rhopalosiphum padi), an important pest of cereal crops, not only directly sucks sap from plants, 
but also transmits a number of plant viruses, collectively the yellow dwarf viruses (YDVs). For quantifying changes in gene 
expression in vector aphids, reverse transcription-quantitative polymerase chain reaction (RT-qPCR) is a touchstone method, 
but the selection and validation of housekeeping genes (HKGs) as reference genes to normalize the expression level of 
endogenous genes of the vector and for exogenous genes of the virus in the aphids is critical to obtaining valid results. 
Such an assessment has not been done, however, for R. padi and YDVs. Here, we tested three algorithms (GeNorm, 
NormFinder and BestKeeper) to assess the suitability of candidate reference genes (EF-la, ACT!, GAPDH, 18S rRNA) in 6 
combinations of YDV and vector aphid morph. EF-1 a and ACTl together or in combination with GAPDH or with GAPDH and 
18S rRNA could confidently be used to normalize virus titre and expression levels of endogenous genes in winged or 
wingless R. padi infected with Barley yellow dwarf virus isolates {BYDV)-PAV and BYDV-GAV. The use of only one reference 
gene, whether the most stably expressed (EF-la) or the least stably expressed (18S rRNA), was not adequate for obtaining 
valid relative expression data from the RT-qPCR. Because of discrepancies among values for changes in relative expression 
obtained using 3 regions of the same gene, different regions of an endogenous aphid gene, including each terminus and 
the middle, should be analyzed at the same time with RT-qPCR. Our results highlight the necessity of choosing the best 
reference genes to obtain valid experimental data and provide several HKGs for relative quantification of virus titre in YDV- 
viruliferous aphids. 
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Introduction 

Rhopalosiphum padi L. (Hemiptera: Aphididae) is an important 
pest of wheat, oat, barley, rye and other gramineous crops in the 
world [1,2]. It can sometimes direcdy affect grain yield by sucking 
nutrients from these plants, but causes the most serious damage 
when it transmits a group of viruses that are collectively called 
yellow dwarf viruses (YDVs) to cereal crops and grasses [3,4]. R. 
padi is considered a competent vector of the YDVs, which 
comprise Barley yellow dwarf virus (BYDV)-PAV and BYDV-PAS 
(genus Luteovirus) and Cereal yellow dwarf virus (CYDV)-RPV, 
CYDV-RPS (genus Polerovirus) and ungrouped BYDV-GPV 
[renamed Wheat yeUow dwarf virus (WYDV)-GPV] [4,5,6,7], aU 
of which belong to the family Luteoviridae [8] . 

R. padi transmits YDVs in a persistendy circulative and 
nonpropagative mode [9,10] after piercing a YDV-infected plant 
and ingesting virions while sucking up the phloem sap [11]. 
Receptor-mediated endocytosis and release must then occur in the 
midgut and/or the hindgut and in the accessory salivary glands 
(ASG) [12,13]. Otherwise, the virions will be excreted along with 



the feces or be retained in the aphid and exposed to a destructive 
environment or immune elements of the aphid or ill-defmed 
metabolic processes so that it cannot complete the circulative route 
in the aphid [14]. The propagation of YDVs in their plant hosts 
has been investigated and elaborated using the prevailing absolute 
and relative quantification methods [15,16], but rarely has 
research focused on fluctuation in the virus titre and accumulation 
in the vector aphid during the acquisition access period (AAP), 
latent period or inoculation access period (lAP), all pivotal issues 
that affect rates of transmission, duration of transmissibility and 
epidemiology of the plant virus in the field. With the entry of the 
virus into the vector aphid, many proteins are incorporated via a 
transcytosis mechanism associated with virus transmission 
[17,18,19,20], followed by up- or down-regulation of the 
expression of numerous genes [21]. Using reverse transcription- 
quantitative polymerase chain reaction (RT-qPCR), we can 
rapidly assess the influence of virus entry on molecular functions 
and biological processes and on the real-time concentration of 
virus [22,2.3]. 
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Table 1. Primers for RT-PCR amplification used in this study. 





Gene 


Forward (F) and reverse (R) primer (5'-3') 


EST size (bp) 


Accession number 


18s rRNA 


F: 


CTGGTTGATCCTGCCAGTAGTCATATG 


571 


KJ612093 


R: TTCCGATTACGGGGCCTCGGATG 


ACT1 


F: 


ATGTGTGACGAAGAAGTAGC 


1131 


KJ612090 


R: TTAGAAGCACTTTCTGTGC 


EF-1a 


F: 


ATTGATATTGCTTTATGGAAATTCG 


837 


KJ612092 


R: ACCAGGGTGGTTCAATACAATAAC 


GAPDH 


F: 


ATGTCAAACATTGGTATCAATGGATTTGG 


999 


KJ612091 


R: TTTAATCCTTAGATTGCATGTACTTGAT 


ago-^a-^ 


F: 


AATCATTTCCAAATTTCAATGCCTCG 


727 


KJ612094 


R: TTGTAACATTACAAACTCTATATTTTC 


ago-la-2 


F: 


ATTAATTCTTTGGTTCGCCGAGCTG 


906 


KJ6 12095 


R: GTACAATATAATTCGATGAGGCTTATAAC 


ago-la-3 


F: 


ATAAAGTTGGAAGGTGATTACAAAC 


278 


KJ612096 


R: TACGTTAAGCATTGTAATTCATCTG 


dcrl 


F: 


TTATGGACTTCCAGACATTAATATTTTATC 


957 


KJ6 12097 


R: TGAGTCATTGCTTGTAATAGATAACTACG 


BYDV-GPV CP 


F: 


ATGAGTACGGTCGCCCTTAGAAATG 


606 




R: CTATTTCGGGTTTTGAAACAGGACC 


BYDV-GPV RTD 


F: 


GTAGACGCGGAACCCGGTCCTAG 


1215 




R: TTAGGAACGGCCACCACCAAATG 


BYDV-PAV CP ' 


F: 


ATGAATTCAGTAGGCCGTAGAGGAC 


600 




R: GGTTCCGGTGTTGAGGAGTCTAC 


BYDV-PAV RTD " 


F: 


GCCAAATAGGTAGACTCCTCAACA 


1341 






R 


: TCACGAAAATTGTGCCTTGTACTC 






R2: TCACTGAAACTGCGCCTTGTATTC 


R3: TCACGAGAATTGAGCCTTGTACGC 


BYDV-GAV CP 


F: 


ATGAATTCAGTAGGCCGTAG 


600 




R: CTATTTGGGAGTCATGTTG 


BYDV-GAV RTD 


F: 


GTAGACTCCTCAACACCAGAG 


1377 




R: TCACCATGTGTCAGCTAAAC 



^With this primer pair, RT-PCR product included less tlian 100 nt of the 5'-terminus of the RTD gene. 
"^With this primer pair, RT-PCR product included less than 100 nt of the 3'-terminus of the CP gene. 
doi:1 0.1 371 /journal.pone.0097038.t001 



A relative quantification method for many aspects of molecular 
research has come into favor to eliminate systemic errors among 
experimental treatments, developmental stages, RNA extraction 
efficiency, and RNA integrity [24,25] . With this method, changes 
in the expression level of genes of interest (GOIs) can be 
determined by comparing the expression of the GOIs with that 
of housekeeping genes (HKGs) [24] . HKGs are constitutive in all 
cells under both normal and disturbed conditions and expressed 
constantly in natural situations [26]. If the expression of HKGs is 
validated as remaining constant under a defined set of conditions, 
relative quantification of GOIs can be realized with convincing 
results [24]. The selection of HKGs for normalization of the 
expression of the GOIs has become routine in studies on human 
disease diagnosis [27,28]. With the pubhcation of the complete 
genomes of different insects, from Drosophila melanogaster [29] to 
Plutella xylostella L. [30], these holometabolous insects have been 
the focus of searches for HKGs to use as reference genes 
[31,32,33]. Estimating the level of expression of HKGs to assay 
differences in expression in hemimetabolous insects, such as 



aphids, which transmit many plant viruses, has started recendy 
with Aphis glycines [34] . For detecting YDV virions in the vector R. 
padi and determining the relative level of expression of vector genes 
using relative RT-qPCR, no HKG has been validated as a suitable 
reference gene to normalize expression. Although ACTl was used 
to determine two YDVs in R. padi, it was actually from a 
coordinate oi Bombyx mori, but the ACTl gene sequence of R. padi 
was unknown [35,36]. 

Selecting appropriate and reliable HKGs as reference genes is 
required for RT-qPCR experiments. Software tools to estimate 
and validate reference genes include GeNorm [25], Normfinder 
[37], and BestKeeper [38]. HKGs such as 18S rRNA, ACTl, EF- 
la, GAPDH, UBI, RPSs, RPLs and TUB are commonly used 
reference genes for standardization when analyzing the expression 
of insect genes using RT-qPCR [34,39,40,41,42,43,44]. 

The announcement of the complete genome of the pea aphid, 
Acyrthosiphon pisum (Harris) [45], has opened a new avenue for 
finding homologous genes in other aphid species. The HKGs oi R. 
padi can be chosen and cloned based on sequences from A. pisum. 
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Figure 1. Box and whisker plots of Cq values for the 4 candidate reference genes in 6 experimental groups (one of three viruses in 
either winged or wingless adults of Rhopalosiphum padh. Each box shows the lower 25* and upper 75* percentiles with median Cq values; 
the whiskers mark the lower 5* and upper 95* percentiles of the Cq values in each data set. Experimental groups from left to right: WYDV-GPV in 
winged adult, BYDV-PAV in winged adult, BYDV-GAV in winged adult, WYDV-GPV and wingless adult, BYDV-PAV and wingless adult, and BYDV-GAV 
and wingless adult. 
doi:l 0.1 371/journal.pone.0097038.g001 



then used as reference genes to quantify the relative titre of viruses 
ingested by R. padi and the relative expression values of 
endogenous genes in the YDV-viruhferous aphid, replacing 
methods for absolute quantification [46], which may lead to 
incorrect, misleading conclusions [15]. 

Previously, actin protein and GAPDH protein of Mjzus penicae 
were found to possess interaction with Beet western yellows virus 
in vitro [18]. In a recent study, some ACT genes (ACTl to ACT4) 
and one GAPDH gene, together with many other genes from A. 
pisum were considered to be potentially related to the transmission 
of Pea enation mosaic virus and Soybean dwarf virus [2 1] . Here, we 
selected these two doubtful HKG genes, ACTl and GAPDH, as 
well as EF-lot and 18S rRNA which were always used as reference 
genes for normalization of target genes in other aphids [42,47], as 
four candidate reference gene for relative RT-qPCR analysis in R. 
padi. Firstly, these four HKGs from R. padi were cloned, and then 
their iiui k'ic a( id scrjuences were aligned with coordinate's from 
other insect species on NCBI. We evaluated the expression 
stability of these four candidate reference genes using 3 tools 
(GeNorm, NormFinder, and BestKeeper) in 6 experimental virus- 
aphid morph combinations after various durations of aphid 
feeding on oat plants infected with the respective virus. To further 
explore the applicability of the optimal reference genes selected by 
GeNorm, we monitored the relative titre of the three viruses in the 
two morphs of R. padi during AAP and compared virus titre 
among the feeding durations using a one-way ANOVA. Also, the 
expression of two endogenous aphid genes, ago-la and dcrl, whose 



encoded proteins are involved in miRNA pathway in aphids and 
have similar functional domains with the RNAi pathway-related 
Ago-2 and Dcr2 proteins [48,49,50], were normalized using the 
selected reference genes to detect any change in expression. 

Results and Discussion 

Sequencing of Candidate Reference Genes from R. padi 

The total RNAs extracted from plants and aphids were all of 
good integrity (Figure SI) and purity (A260/A280 values: 1.95 ~ 
2.10; A260/A230 values > 2.0). All genes were amplified by RT- 
PCR (Figure S2) with the primers described in Table 1. Nucleic 
acid sequences of 4 candidate reference genes yielded good 
matches in a BLAST search of the NCBI database with high 
similarity to corresponding genes from other aphids and insects: 
18S rRNA, 90-100%; EF-la, 83-94%; ACTl, 82-96%; 
GAPDH, 78-95%, respectively (Table SI). 

Expression Profiles of Candidate Reference Genes 

In the two-step RT-qPCR, the amplification efficiencies of the 4 
candidate reference genes ranged from 91% to 97% (Table 2; 
Figure S3). Relative standard cur\'es for determining the 
amplification efficiencies of these genes are shown in Figure S4, 
and all Cq values of the dilution series of the cDNA for each gene 
were less than 30, while those of the no-template controls were 
either higher than 36 or undetermined. Consequentiy, the raw Cq 
value of the individual gene reflected the actual mRNA level of this 
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Table 2. Primers for RT-qPCR amplification used in this study. 





Gene 


Forward and reverse primer sequence* 


Amplicon 


Efficiency' 






size" (bp) 


(winged/wingless) 


18S rRNA 


5'-CGGGAGGAACGCTTTTATTAGA-3' 


262 


90.506/91.014 




5 ' -TTGG ATGTGGTAGCCGTTTCTC-3 ' 






ACTl 


5'-AACGGAAGCACCTTTGAACC-3' 


385 


93.665/95.067 


5'-GGAAGAAGCAGCAGTAGCCAT-3' 


EF-lo! 


5 ' -TAG ACGCTATCCTACCCCCC A-3 ' 


328 


95.902/97.157 


5 ' -GTG A A ATC AGC AGC ACCCTTG-3 ' 


GAPDH 


5 '-ACTACTGTTCACGCTACCACCG-3 ' 


272 


94.814/96.216 


5 ' -GCTGCTTCCTTG ACCTTACCTT-3 ' 


ago-!a-l 


5 ■ - ACC AGTCTGTTCGTCC ATCTC AAT-3 ' 


138 


91.891/91.594 


5 ' -GTTTTCTTTGTTC ACC A ATGTCCC-3 ' 


ago-la-2 


5 '-ACCG AGCATCAGACCTAAAGTATT-3 ' 


215 


98.391/99.05 


5'- TTAGAAGTTCCCTGACCATAGAGC-3' 


ago-la-3 


5'-TATTCTGTGCCGACAAAAAGG-3 ' 


184 


97.263/91.749 


5 ' -G ATTC A A AATGGTTGTC ATCCC-3 ' 


dal 


5'-CTACCACCCTGTTACTTTGTCCC-3' 


168 


93.875/101.698 


5'-GTCGTCATTTGAGAAGCCCAC-3' 


BYDV-GPV CP 


5'-TTCGTTTTCGCAAAGGATTCA-3' 


359 


94.75/101.875 


5'-GTG"n"GGCGTCACCG"n"ACC-3' 


BYDV-GPV 


5 ' -G ACCATC ACCCACTCC ACCT-3 ' 


369 


102.875/96.203 


RTD 


5 ' -CCTCGCTTG A ATC ATC ATACG-3 ' 






BYDV-PAV CP 


5 '-CGGGGCTGAGGTATTCGTAT-3 ' 


281 


93.65/98.105 


5 ' - AGG ACTTTG AGGCGG ATTTG-3 ' 


BYDV-PAV 


5 ■ -C ATTGCCTACTCC AACTC ATCG-3 ' 


288 


100.303/92.586 


RTD 


5 ' -ATTCTTTTGTTCGTGTACCCTCC-3 ' 






BYDV-GAV CP 


5 '-CAGGCAGG ACTG AGGTATTCGTA-3 ' 


278 


96.752/100.563 


5'-GGTTGCTGA"rnTGAGATGGTGA-3' 


BYDV-GAV 


5 ■ -GTGCTG ACCCC A AGTTTTACCTC-3 ' 


191 


98.585/100.127 


RTD 


5 ' -GTTTTCTGTTTTCCTAACCGTGCT-3 ' 







"Designed by an automated search in Primer Premier 5 software. Parameters for each primer pair: C+C content = 40-60%, annealing temperature between 58° and 
62"C, and no false priming. 

"^Amplicon size and specificity were checked using melt curve and 2% agarose gel electrophoresis (Figure S3). 
'Efficiency % = 1o'""''°f"" -1. fi^ values of all relative standard curves were more than 0.989. 
doi:1 0.1 371 /journal.pone.0097038.t002 



gene and could be compared directly among the 6 experimental 
groups. The distribution of the Cq values of the 4 candidate 
reference genes showed that the expression of 18S rRNA was 
maintained at extremely high levels with the lowest Cq values 
(Cq= 12~14), higher levels with moderate Cq values (Cq 16~22) 
were obtained for ACTl, EF-la and GAPDH, but EF-la was at 



the lower edge of this range in all experimental groups (Figure 1). 
The high expression level of 18S rRNA in R. padi was expected, in 
agreement with the lowest Cq value in other insects [44] . EF- 1 ot, a 
weU-known HKG in various organisms [51] and an evolutionarUy 
conserved nuclear gene in aphids [52,53,54], can be used to 
classify the population oi R. padi into various subgroups depending 



Table 3. Influence of YDV entry and wing morph on raw Cq value for individual genes. 



Gene Virus Wing morph interaction 



18S rRNA 


0.000*** 


0.285"' 


0.276"' 


ACTl 


0.000*** 


0.040* 


0.000*** 


EF-la 


0.000*** 


0.874"' 


0.504"' 


GAPDH 


0.010* 


0.000*** 


0.020* 



Two-way ANOVA analysis: ns, nonsignificant, p>0.05; * p<0.05, ** p<0.01; *** p<0.001. 
doi:10.1371/journal.pone.0097038.t003 
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Table 4. Paired connparisons for reference genes between the WYDV-GPV (GPV), BYDV-PAV (PAV) and BYDV-GAV (GAV)- 
viruliferous aphids {Rhopalosiphum padi). 



Gene 


GAV 




GPV 




PAV 




GPV 


PAV 


GAV 


PAV 


GAV 


GPV 


18S rRNA 


0.000*** 


0.000*** 


0.000*** 


0.879 


0.000*** 


0.879 


ACT1 


0.000*** 


0.310 


0.000*** 


0.000** 


0.310 


0.000*** 


EF-la 


0.000*** 


0.139 


0.000*** 


0.000** 


0.139 


0.000*** 


GAPDH 


0.009** 


0.590 


0.009** 


0.109 


0.590 


0.109 



Tukey's HSD post hoc tests: ns, nonsignificant, p>0.05; * p<0.05; ** p<0.01; *** p<0.001. 
doi:1 0.1 371 /journal.pone.0097038.t004 



on EF-la diversity in the population [52,53]. Its expression level 
was consistent between winged aphids and wingless aphids, lower 
than that of 18S rRNA but slighdy higher than that of ACT 1 and 
GAPDH. The expression levels of ACTl and GAPDH were 
similar in R. padi, but the mRNA level of ACTl largely differs 
from GAPDH in planthoppers [44]. ACTl and GAPDH encode a 
cytoskeletal structural protein and a catalytic enzyme, respectively 
[44]. Thus, it was difficult to interpret here why these two genes 
with diverse functions actually have similar expression levels in R. 
padi. 

The influence of virus and aphid morph on the Cq values of the 
4 candidate reference genes was statistically analyzed with a two- 
way ANOVA (Table 3). When comparing Cq values for each gene 
(18S rRNA, ACTl and EF-la) among the different YDV- 
viruliferous aphids, Cq differed significantly (all /)<0.001); values 
for GAPDH differed less among the three viruses (0.01</)<0.05). 
For the aphid morph factor, wing morph had no significant effect 



on the expression of 18S rRNA (/; = 0.285) or EF-la (p = 0.874), 
but did contribute to the significant difference in ACTl (/) = 0.04) 
and in GAPDH (/)<0.001). For the interaction of YDV and wing 
morph, Cq values of 18S rRNA and EF-la revealed no significant 
difference (both p>0.05), as opposed to the Cq values for ACTl 
(/)<0.001) and GAPDH ip = 0.02), probably due to the significant 
differences for both virus and wing morph. 

Since the virus was considered the main effect and we had more 
than 2 viruses, we used Tukey's honestly significant difference 
(HSD) post-hoc test to compare the influence of (1) BYDV-GAV 
vs. WYDV-GPV, (2) BYDV-GAV vs. BYDV-PAV, and (3) 
WYDV-GPV vs. BYDV-PAV on each candidate reference gene 
(Table 4). For the influence of infection by BYDV-GAV vs. 
WYDV-GPV on the expression of 18S rRNA, ACTl or EF-la of 
R. padi, the expression levels of these candidate reference genes 
differed significantly (all p<0.00l), and less so for GAPDH 
(/) = 0.009). Differences in 18S rRNA expression for BYDV-GAV 



Table 5. Stability ranl<ings according to three software tools for individual endogenous reference genes in winged adults of 
Rhopalosiphum padi carrying one of three YDVs. 





Virus and reference gene 


GeNorm * 


NormFinder'' 


BestKeeper*^ 




M 


stability value 


Cq set SD 


GPV 


EF-la 


0.396 


0.080 


0.29 


ACTl 


0.406 


0.085 


0.34 


GAPDH 


0.421 


0.110 


0.30 


18S rRNA 


0.629 


0.192 


0.45 


PAV 


ACTl 


0.270 


0.054 


0.31 


EF-la 


0.274 


0.049 


0.31 


GAPDH 


0.309 


0.094 


0.28 


18S rRNA 


0.444 


0.165 


0.26 


GAV 


EF-la 


0.211 


0.051 


0.29 


ACTl 


0.226 


0.074 


0.31 


GAPDH 


0.259 


0.076 


0.27 


18S rRNA 


0.338 


0.151 


0.21 



^According to GeNorm, the smallest M value indicates the most stable gene expression, the largest M value the most unstable; expression of a gene with M >1.5 is 
considered inconsistent. 

"^Like GeNorm, the most stable gene has the lowest stability value, and vice versa. 

■^Genes with SD <1 can be ranked according to their Cq SD values. The gene most stably expressed has the lowest SD value; the highest SD value indicates the gene 
with the most unstable expression. 
doi:1 0.1 371 /journal.pone.0097038.t005 
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Table 6. Stability rankings according to three software tools for individual endogenous reference genes in wingless adults of 
Rhopalosiphum padi carrying one of three YDVs. 





Virus and reference gene 


GeNorm* 


NormFinder'' 


BestKeeper*^ 




M 


stability value 


Cq set SD 


GPV 


EF-la 


0.508 


0.033 


0.26 


GAPDH 


0.570 


0.115 


0.46 


1 8S rRNA 


0.648 


0.202 


0.29 


ACT1 


0.799 


0.231 


0.68 


PAV 


EF-la 


0.277 


0.017 


0.36 


ACT1 


0.285 


0.028 


0.38 


GAPDH 


0.296 


0.060 


0.34 


18S rRNA 


0.535 


0.213 


0.24 


GAV 


EF-la 


0.293 


0.053 


0.42 


ACTl 


0.294 


0.054 


0.49 


GAPDH 


0.317 


0.078 


0.48 


18S rRNA 


0.527 


0.212 


0.29 



^According to GeNorm, the smallest M value indicates the most stable gene expression, the largest M value the most unstable; expression of a gene with M >1.5 is 
considered inconsistent. 

"^Like GeNorm, the most stable gene has the lowest stability value, and vice versa. 

'^Genes with SD <1 can be ranked according to their Cq SD values. The gene most stably expressed has the lowest SD value; the highest SD value indicates the gene 
with the most unstable expression. 
doi:l 0.1 371 /journal.pone.0097038.t006 



VS. BYDV-PAV was highly significant (/)<0.001), but not 
significant ior the other 3 reference genes (ACTl, /> = 0.310; EF- 
la, j6=0.139; GAPDH,/) = 0.590). WYDV-GPVand BYDV-PAV 
difiered significantly in expression of ACTl and of EF-la (both 
^<0.001), but consistendy lacked difierences for 18S rRNA 
(/! = 0.879) and for GAPDH (/;= 0.109). 

In summary, the Cq values for expression of 18S rRNA in 
WYDV-GPV- and BYDV-PAV-viruliferous R. padi fell into one 
group (subset 2), the BYDV-GAV-viruliferous R. padi into another 
(subset 1). As for the Cq values of ACTl and EF-la, die WYDV- 
GPV-viruliferous R. padi was in one group (subset 1), the others in 
a second (subset 2). The Cq values for die GAPDH in BYDV- 
PAV- and BYDV-GAV-viruliferous R. padi are logically within 
one group, rather than with WYDV-GPV-viruliferous R. padi, 
because of the lack of significant differences between the infection 
with BYDV-PAV and BYDV-GAV (/; = 0.590), with /)-values 
lower than that for the comparison between WYDV-GPV and 
BYDV-PAV (/)= 0.109). The lack of unanimous grouping of Cq 
values for these 4 candidate reference genes might be attributed to 
the differing transmission properties of the YDVs by R. padi 
because successful transmission depends on the compatibility 
between the virus and vector insects [9] . 

Expression Stability Analysis of Reference Genes in YDV- 
viruliferous R. padi 

For analyzing the expression stability of the 4 candidate 
reference genes, a stability value for each individual gene was 
calculated using three Excel-based tools (GeNorm [25] and 
NormFinder [37] and BestKeeper [38]). The amplification 
efficiencies of 4 candidate genes in winged R. padi differed 
somewhat from that of the wingless (Table 2). Therefore, the 
stability ranking fell into two sets. 



GeNorm only recognizes linear scale quantities, obtained by 
transforming raw Cq values using the delta-Cq method (described 
in the methods) or a standard curve; an expression stability value 
M for all genes contained in the InputData file is automatic 
calculated by an algorithm [25] and displayed. The smallest M 
value indicates the most stable gene expression; the largest Af value 
indicates the most unstable expression. Importantiy, expression of 
a gene with Af > 1 .5 is considered to be inconsistent [25] . After the 
Cq values for the RT-qPCR of all 4 candidate reference genes 
under the 6 experimental conditions were transformed into a valid 
data set (see formula in the methods), the M values were calculated 
by a GeNorm applet. As shown in Tables 5 and 6, 4 virus-aphid 
combinations (WYDV-GPV-winged, BYDV-GAV-winged, 
BYDV-PAV-wingless and BYDV-GAV-wingless) had the same 
order of expression stability from the most stable to the least stable 
gene: EF-la, ACTl, GAPDH and 18S rRNA. For the BYDV- 
PAV-winged combination, ACTl was ranked as the best reference 
gene (lowest M), followed by EF-la, GAPDH and 18S rRNA, 
However, ACTl was the worst in the case of WYDV-GPV- 
wingless, and the stability sequence of the other 3 genes was the 
same for BYD V-PAV-winged. In brief, EF- 1 a (5 first places and 1 
second place) seemed to be the best reference gene to normalize 
the GOIs in aU YDV-viruliferous types of R. padi. 

A different algorithm, NormFinder, only analyzes linear scale 
data [37]. The stability values of the reference genes can be 
calculated for one experimental group or for multiple groups in 
one data set. The most stable gene has the lowest stability value, 
and vice versa [37]. When the transformed linear data (the same Q_ 
values of the 4 candidate reference gene analyzed by GeNorm) 
was analyzed with NormFinder, sorting for expression stability of 
the candidate reference genes (Tables 5 and 6) for the WYDV- 
GPV-winged, BYDV-PAV-winged, BYDV-GAV-winged and 
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Figure 2. Determination of optimal reference genes for normalization by GeNorm. The optimal number of reference genes was 
determined by pairwise variation V (A) between two sequential normalization factors containing an increasing number of reference genes (n and n +^ 
reference genes). The pairwise variation begins with the first two and three (1/2/3) most stable candidate reference genes (B), terminating with the 
addition of the least stable one, here 1/3/4. The cut-off 1/ value is set at 0.15 by the program GeNorm, below which an additional reference gene does 
not need to be included for calculating the normalization factor. 
doi:1 0.1 371 /journal.pone.0097038.g002 



BYDV-PAV-wiiigles.s groups was highly uniform: EF-la was the 
most stable gene, but 18S rRNA was the most unstable. EF-lot was 
also the first choice of reference genes for the WYDV-GPV- 
wingless group, but sorted as the third most stable gene for the 
BYDV-GAV-wingless. 



Finally, the raw Cq values of candidate reference genes in each 
group were assessed directly by BestKeeper; when the SD of the 
Cq is higher than 1, expression for that is considered as 
inconsistent, and only genes with a SD less than 1 can be ordered 
using their Cq SD values [38]. The gene most stably expressed 
thus has the lowest SD value; the gene most unstably expressed has 
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Figure 3. Mean relative expression values (± SE, n=3) for endogenous aphid genes ago-7a-1 (A), ago-1a-2 (B), ago-7a-3 (C) and dcr7 
(D) in BYDV-GAV-viruliferous winged adults of Rhopa/osiphum padi after different virus-feeding durations. A-D, Expression values for 
each gene or each region of one gene at each duration were normalized with the reference gene(s) selected by GeNorm and then compared with a 
one-way ANOVA (p) among these durations with each normalization condition. 2re = 2 best reference genes; 3re = 3 best reference genes; 4re = all 4 
reference genes; bad2re = 2 least stable reference genes; EF-1a = EF-1o( as the reference gene; 18S = 18S rRNA as the reference gene. 
dol:1 0.1 371 /journal.pone.0097038.g003 



the highest SD value. In this study, Cq SD values from the qPCRs 
of the 4 candidate reference genes in all groups were less than 1 
(Tables 5 and 6) and can be compared with each other as 
individual groups. Unlike the results from GeNorm and 
NormFinder, 18S rRNA was evaluated as the first and second 
most stable gene in the following 5 groups: BYDV-PAV-winged, 
BYDV-GAV-winged, BYDV-PAV-wingless, BYDV-GAV-wing- 
less and WYDV-GPV-wingless, but as the least stably expressed 
gene in the GPV-winged group. EF-lot, assessed as a nearly ideal 
reference gene by GeNorm and NormFinder tools, was ranked as 
the optimal reference gene only in the WYDV-GPV-winged and - 
wingless groups. For the other 4 groups, the performance of EF-la 
was mediocre. 

As discussed already, GeNorm and NormFinder yielded very 
similar results with respect to sorting expression stability in each 
group and in identifying the best reference gene among 4 
candidate reference genes. We inferred that the two programs 



analyzed the same linear scale data set, which was transformed 
from raw Cq values and corrected by amplification efficiency of 
the respective gene (see formula in the methods). The results from 
the BestKeeper tool were nearly opposite to those obtained with 
GeNorm and NormFinder. Perhaps the number and variety of 
candidate reference genes that we used here are not very 
diversified. 

Selection of Reference Genes for Normalization 

There were 3 cases for calculating the best number of reference 
genes in our 6 experimental groups: (1) F2/3<0.15, F3/4<0.15; 
(2) F2/3<0.15, F3/4>0.15; and (3) F2/3>0.15, I'3/4>0.15 
(Figure 2A). The first case meant the best number of reference 
genes could be 2, 3 or 4; simply, the 2 most stable genes, ACTl 
and EF- 1 a, (Figure 2B) were the ideal combination of reference 
genes for the normalization analysis of target genes in the 
respective groups, BYDV-PAV-winged, BYDV-GAV-winged, 
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BYDV-PAV-wingless, and BYDV-GAV-wingless. The second 
case demonstrated that the optimal choice of reference genes 
was ACTl and GAPDH for the WYDV-GPV-winged group 
(Figure 2B), while adding another one or two candidate reference 
genes increased the risk of misinterpreting the expression profile of 
the GOIs. The last case indicated that using just thes(; 4 candidate 
reference genes would not form a dependable reference gene 
group for the relative quantification of GOIs in WYDV-GPV- 
viruliferous wingless R. padi. 

Relative Quantification of Endogenous Genes 

Since no up- or down-regulation of these candidate internal 
control genes was observed in the 6 combinations (SD <1 from 
BestKeeper and stability value M <1.5 from GeNorm), and the 
optimal number of reference genes had been selected by GeNorm, 
these genes can be used in the expression profile analysis of 
endogenous and low-abundance genes. Because the expression 
levels of ago-la and dcrla (dcrl) gene varied significandy among 
four reproductive morphs of ^. pisum when measured with semi- 
quantitative RT-PCR [49], miRNA machinery and alternative 
transcription are thought to play a key role in morph transforma- 
tion in the life cycle of A. pisum. 

Regardless of the normalization method used (described in 
Materials and Methods), no significant difiFerences were observed 
in the expression level for the 3 regions of the ago-la and dcrl genes 
before and after the entry of BYDV-GAV into winged (Figure 3 
A-D) or wingless (Figure 4 A-D) adults of R. padi (all /)>0.05, 
except for when EF-loi was used as the reference gene to 
normalize the expression change of ago-la-3 region \p — 0.038] 
[Figure 3C]). We thus concluded that the expression of the ago-la 
and dcrl genes did not respond to the entry of BYDV-GAV into R. 
padi and remained unchanged. This finding conformed to the facts 
that BYDV-GAV could scarcely be attached by receptor proteins 
of R. padi, and has rare chance to access to hindgut apical 
plasmalemma [12]. 

The entry of BYDV-PAV did not seem to affect the expression 
profile of ago-la in either of the wing morphs oiR. padi on the basis 
of the results in Figures 5A and 6A, for relative expression values 
of ago-la-l among these feeding durations approached the same 
(all ^>0.05 under all normalizations). However, for expression 
levels o{ ago-la-2 and ago-la-3 (the middle region and 3'-terminus 
of ago-la gene [49]) using all the same normalization methods, we 
drew difierent conclusions: after being normalized with a single 
reference gene (EF-la or 18S rRNA) or with bad reference genes, 
the expression of ago-la-2 and -3 in BYDV-PAV-viruliferous 
winged R. padi showed no significant changes from 0 h to 96 h 
feeding duration (purple, yellow and red bars in Figure 5B and C, 
all ^>0.05); thus, ago-la-2 and -3 appeared to be unaffected by 
BYDV-PAV acquisition; when the best 2 or 3 reference genes or 
all 4 reference genes were used to normalize the expression of 
these two regions in winged R. padi, a difference in the expression 
of ago-la-2 among the diflFerent feeding durations emerged slightiy 
above the significance level (blue, green, and brown bars in 
Figure 5B, 0. 1>^S0.05), but the expression values for ago-la-3 
already declined significandy before (0 h) and immediately after 
entry of BYDV-PAV (blue, green, and brown bars in Figure 5C, 
^<0.().5) in sj)itc of recover to initial level afterwards. In wingless R. 
padi, the changes in expression for ago-la-3 higher than that for ago- 
la-2 (Figure 6B and C, both/)<0.1, the p value for ago-la-3 was 
lower than for ago-la-2) was observed when the best 2 or 3 
reference genes or all 4 reference genes were used as normalization 
methods, with a similar trend in the expression change for ago-la-2 
and ago-la-3 in the winged morph. On the basis of these findings, 
we concluded that although no change in expression was detected 



for ago-la-1 in R. padi, the expression of the full-length ago-la was 
disturbed by the entry of BYDV-PAV particles as seen by the 
obvious expression change in ago-la-2 and ago-la-3. For the 

expression change in dcrl in BYDV-PAV-viruliferous winged and 
wingless R. padi (Figures 5D and 6D), the significance level (both 
p<0.l) was near that of ago-la-2 and -3 (all p<0.l) with the 
normalization factors calculated by the best 2, or 3 reference genes 
or all 4 reference genes. The detected region for the dcrl gene is 
also in the middle region and the 3 '-terminus of the full-length dcrl 
gene [49], so we speculated that the entry of BYDV-PAV also 
interfered with the expression oi dcrl. The unexpected discrepan- 
cies among values for changes in relative expression obtained using 
3 regions of the same gene in BYDV-PAV-viruliferous R. padi 
suggested that mechanisms on the splicing and reconstruction of 
these three regions in ago-la were difierences in aphids, probably 
because the pre-mRNA of region 2 and region 3 were recruited 
and then translate into the PAZ and PIWI domains of antivirus 
defense-associated AG02 protein in aphids. 

For WYDV-GPV-viruliferous R. padi, we could not give an 
explicit inference that the relative expression in the 3 regions of 
ago-la and dcrl gene changed or unchanged with the virus 
acquisition, since these normalization methods provided no 
concurrent and regular results (Figures 7 and 8). We inferred that 
the complex interactions between WYDV-GPV, oat plant and 
vector R. padi led to the unsuitability of the selected reference genes 
(Figure 2 A, V2/3 and F3/4 near or above the cut-off value of 
0.15). 

The miRNA pathway has close relationship with RNAi 
pathway in eukaryotes, and the two pathways possess extremely 
parallel mechanism of RNA cleavage in cell cytoplasm [.^5] . The 
machinery proteins of miRNA pathway. Ago- land Dcrl, have 
high similarity with those of RNAi pathway, Ago-2 and Dcr2 
[49,56]. Since R. padi is eHicient vector of BYD\'-P.W and 
WYDV-GPV, these RNA pathways may be triggered against virus 
invasion [56]. 

Fluctuation and Comparison of Virus litre in /?. padi 

Fluctuation in the virus titre in R. padi was monitored by two 
methods. The first method used the SD values of the raw Cq 
values (Cq SD) for the CP gene and RTP gene (the Cq values of 
the no-template controls and 0-h-feeding duration sample of these 
two genes were either higher than 36 or undetermined and thus 
not analyzed here) as calculated by BestKeeper. In WYDV-GPV 
and BYDV-PAV, viruliferous winged and wingless R. padi, the Cq 
SD values of the CP gene and RTD gene were all higher than 1, 
which was affirmed to be inconstant and statistically significandy 
different among feeding durations; but the Cq SD values of the CP 
gene and RTD gene of BYDV-GAV in winged and wingless R. 
padi were lower than 1, with no variance in virus titre during the 4- 
day AAP (Table S2 and S3). These results agreed with the 
transmission specificity of these three YDVs by R. padi, and the 
difference in Cq at the initial virus-acquisition duration (maximal 
Cq value) implied the presence of a sluggish receptor-mediated 
recognition process at the midgut and/or hindgut [9]. One 
problem with this method is that the trend for increased virus load 
is inversely related to the duration of virus-feeding. 

Tlu' second one is the relative virus load, normalized l)y one 
reference gene or a combination of referencx genes (with the same 
normalization as the endogenous genes). The expression levels of 
tiie CP gene and RTD gene of WYDV-GPV, BYDV-PAV and 
BYDV-GAV in winged (Figure 9 A-F) and wingless R. padi 
(Figure 10 A-F) increased more than 5-fold from the extraordi- 
narily low value at the initial 1 2 h to the highest value during the 
last day (72 to 96 h), and the difference in expression of either the 
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Figure 4. Mean relative expression values (± SE, n=3) for endogenous aphid genes ago-7a-1 (A), ago-7a-2 (B), ago-la-3 (C) and dcr7 
(D) in BYDV-GAV-viruliferous wingless adults of Rhopa/osiphum padi after different virus-feeding durations. A-D, Expression values for 
each gene or each region of one gene at each duration were normalized with the reference gene(s) selected by GeNorm and then compared with a 
one-way ANOVA (p) among these durations with each normalization condition. 2re = 2 best reference genes; 3re = 3 best reference genes; 4re = all 4 
reference genes; bad2re = 2 least stable reference genes; EF-loi = EF-lo( as the reference gene; 18S = 18S rRNA as the reference gene. 
doi:l 0.1 371 /journal.pone.0097038.g004 



CP or RTD gene across all feeding durations was significant (all 
p<0.05, Figures 9 A-F and 10 A-F), regardless of the 
normalization setting adopted. When tlie duration of feeding on 
either WYDV-GPV-, BYDV-PAV- or BYDV-GAV-infected oat 
plants exceeded 48 h, the relative expression values of the CP gene 
and RTD gene continued to rise at a slow rate and rose and fell 
within a narrow range. The findings reflected accumulation but no 
tendency for replication of the YDV particles in R. padi and the 
saturability of virus load by aphid, in agreement with Gray and 
Gildow's hypothesis about the nonpropagative transmission of 
YDVs by aphids [9] . As shown in Figures 9 and 1 0, the fluctuation 
in virus litre of the three viruses in winged and wingless R. padi at 
each duration within the 4-day AAP had some slight differences, 
relied on the method of normalization, and a good reference gene 
or reference gene combination estimated by GeNorm resulted in a 
much lower mean SE value for the sample replicates. 



We could not uncover differences in transmissibility of these 
three YDVs by R. padi by comparing the relative virus tilre during 
the first 2 days of feeding with a one-way ANOVA because the 
relative virus litre between transmissible WYDV-PAV or BYDV- 
PAV and non- transmissible BYDV-GAV did not differ signifi- 
cantly) (all p>0.05, data not shown), even though the relative titre 
of WYDV-GPV and BYDV-PAV (Figures 9 A-D and 10 A-D) for 
the 12, 24 and 36 h of feeding was lower than that of BYDV-GAV 
(Figures 9E, 9F, lOE, and lOF). Perhaps the typical biological 
transmission assay is the most efficient and convincing measure for 
determining the insect vector's transmission specificity for a 
particular plant virus [19]. 

With the extension of aphid feeding duration on infected oat 
plants, the relative expression values of the CP or RTD gene of all 
three YDVs in R. padi increased gradually from 1 2 h to 60 h or 
72 h, but tended to stabilized after 60 h or 72 h. BYDV-PAV and 
WYDV-GPV once acquired are retained life of R. padi. Following 
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Figure 5. Mean relative expression values (± SE, n=3) for endogenous aphid genes ago-7a-1 (A), ago-7a-2 (B), ago-1a-2 (C) and dcrl 
(D) in BYDV-PAV-viruliferous winged adults of Rhopalosiphum padi after different virus-feeding durations. A-D, Expression values for 
each gene or each region of one gene at each duration were normalized with the reference gene(s) selected by GeNorm and then compared with a 
one-way ANOVA (p) among these durations with each normalization condition. 2re = 2 best reference genes; 3re = 3 best reference genes; 4re = all 4 
reference genes; bad2re = 2 least stable reference genes; EF-1a = EF-1o( as the reference gene; 18S = 18S rRNA as the reference gene. 
doi:1 0.1 371 /journal.pone.0097038.g005 



ingestion during feeding, the virus moves through the alimentary 
canal of R. padi, invades the gut epithehal cells and then is released 
into the hemolymph, ultimately infects the salivary glands and is 
released into the salivary ducts where it can be transferred to new 
plants via the saliva released during feeding [57]. Virus litre 
continues to rise during the first 48 h of the AAP and lAP, but 
decreases slightly because the virus particles are transmitted into 
plants. As aphids continue sucking the virus-containing phloem 
sap, the virus load in aphids rises again. For the untransmissible 
one, the virus only moves through the alimentary canal and then 
direct excretes from body. Therefore, the virus titre increases at 
initial feeding duration, and then keeps a dynamic balance because 
of the continuous sucking [12,58]. 

Conclusions 

In this study, we found that in testing EF-lot, ACTl, GAPDH 
and IBS rRNA as reference genes, combining the first 2 or 3 genes 



or using all 4 genes can be used to normalize the relative 
expression level of endogenous and virus titre in BYDV-PAV and 
BYDV-GAV-viruliferous winged and wingless R. padi. Using only 
one reference gene, regardless of its expression stability, is not 
recommended for relative quantification by RT-qPCR. Impor- 
tantly, for accurate and precise results, the use of different regions 
of a gene, including the 5'-terminus, middle and 3'-terminus, at 
the same in the relative RT-qPCR is strongly advised, since the 
different regions of one gene endure different selection pressures 
and may evolve some other new functions in complex life cycle of 
aphids. To assess relative expression and virus load in winged and 
wingless WYDV-GPV-viruliferous R. padi, the optimal number of 
combination of reference genes should be selected carefully, 
abandoning genes that are not useful alone or in combination and 
testing new ones as needed. 
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Figure 6. Mean relative expression values (± SE, n = 3) for endogenous aphid genes ago-1a-^ (A), ago-1a-2 (B), ago-1a-2 (C) and dcrl 
(D) in BYDV-PAV-viruliferous wingless adults of Rhopalosiphum padi after different virus-feeding durations. A-D, Expression values for 
each gene or each region of one gene at each duration were normalized with the reference gene(s) selected by GeNorm and then compared with a 
one-way ANOVA (p) among these durations with each normalization condition. 2re = 2 best reference genes; 3re = 3 best reference genes; 4re = all 4 
reference genes; bad2re = 2 least stable reference genes; EF-1oi = EF-1o( as the reference gene; 18S = 18S rRNA as the reference gene. 
doi:1 0.1 371 /journal.pone.0097038.g006 



JVIaterlals and IVIethods 

Virus Maintenance and Aphid Rearing 

Laboratory isolates of XWDV-GPV, BYDV-PAV, and BYDV- 
GAV have been maintained on oat plants [Avena sativa K. Koch cv. 
Coast-Black Oat) in our laboratory since the 1990s [59]. A virus- 
free laboratory population of R. padi was reared on winter wheat 
(Triticum aestivum L., cv. Yangmai 158), under controlled conditions 
at 18-23°C. A low-density population was maintained for the 
wingless morph of R. padi, and a high-density population was 
maintained to obtain the winged morph. This R. padi population 
transmits BYDV-PAV (70-80% transmission rate) and WYDV- 
GPV (50% transmission rate) effectively, but poorly transmits 
BYDV-GAV (5-7% transmission rate) [60]. 

To confirm the morph features of wingless aphids and winged 
aphids, we observed different developmental stages o{ R. padi from 
frrst to fifth instar daily with a stereomicroscope (Olympus 



SZXI6). Because the lifespan of the winged and wingless adult 
aphids (more than 5 days) is much longer than that of each nymph 
instar (less than 2 days) and the adult aphids are distinguishable 
from the nymphs, winged and wingless adults were used for the 
experiments. 

Experimental Treatment and Sampling 

Approximately 300 wingless and 300 winged adult aphids were 
picked from wheat with a soft brush, and then released carefully 
onto two oat plants infected with WYDV-GPV, BYDV-PAV or 
BYDV-GAV for acquisition of viruses, with insect-proof nets 
separating the three sets of infected plants. After each feeding 
duration (0, 12, 24, 36, 48, 60, 72, 84 or 96 h), 24 wingless adult 
aphids in each net were collected and divided randomly into 3 
pools, then 8 aphids of each type were placed into an RNase-free 
tube. The winged adult aphids were treated the same way. The 
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Figure 7. Mean relative expression values (± SE, n = 3) for endogenous aphid genes ago-la-^ (A), ago-Ja-2 (B), ago-la-Z (C) and dcrl 
(D) in BYDV-GPV-viruliferous winged adults of Rhopalosiphum padi after different virus-feeding durations. A-D, Expression values for 
each gene or each region of one gene at each duration were normalized with the reference gene(s) selected by GeNorm and then compared with a 
one-way ANOVA (p) among these durations with each normalization condition. 2re = 2 best reference genes; 3re = 3 best reference genes; 4re = all 4 
reference genes; bad2re = 2 least stable reference genes; EF-1a = EF-1o( as the reference gene; 18S = 18S rRNA as the reference gene. 
doi:1 0.1 371 /journal.pone.0097038.g007 



tubes were immediately placed in liquid nitrogen and stored in a — 
70°C freezer for RNA extraction. 

Selection of Candidate Reference Genes and Target 
Genes and Design of Primers 

Four HKGs ofR. padi, including 18S rRNA, ACTl, EF-lot, and 
GAPDH were selected as candidate reference genes (Table 7). 
Based on the corresponding sequences of R. padi and A. pisum from 
the NCBI GenBank, primers for reverse transcription-PC R (RT- 
PCR) amplification of 18S rRNA and EF-la were designed 
(Table 1). Primer pairs for amplifying the ACTl and GAPDH 
sequences oi R. padi were designed using homologous sequences of 
A. pisum deposited in GenBank (Table 1). 

Two endogenous genes oi R. padi, ago-la and dcrl, which might 
be involved in the miRNA pathway in aphids [48,49], were used 
as target genes. The 3 packed regions of the ago-la have been 
showed different selective pressures, region 1 and region 2 



reflecting strong purifying selection (negative selection), but region 
3 owning no positively selected sites [49]. The chosen region of 
dcrl has the same selection pressure with region 3 of ago-la [49]. 
For the RT-PCR of these two genes, primers (Table 1) were 
devised using the DNA regions for the ago-la gene and dcrl gene. 

The CP gene (ORF3), encoding the coat protein, and RTD 
gene (ORF5), encoding the readthrough protein, of WYDV-GPV, 
BYDV-PAV and BYDV-GAV from infected oat plants were 
amplified by RT-PCR using their respective primers (Table 1). All 
primers used in this study were synthesized by Sangon Biotech 
(Shanghai) Co. 

RNA Extraction and Purification 

For plant leaves and a large number of aphids, samples were 
first homogenized in liquid nitrogen with a sterile mortar and 
pestle. Then total RNA was extracted from the resulting fine 
powder using the instructions for TRIzol reagent (Ambion, USA) 
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Figure 8. Mean relative expression values (± SE, n = 3) for endogenous aphid genes ago-la-^ (A), ago-Ja-2 (B), ago-la-Z (C) and dcrl 
(D) in BYDV-GPV-viruliferous wingless adults of Rhopalosiphum padiaUer different virus-feeding durations. A-D, Expression values for 
each gene or each region of one gene at each duration were normalized with the reference gene(s) selected by GeNorm and then compared with a 
one-way ANOVA (p) among these durations with each normalization condition. 2re = 2 best reference genes; 3re = 3 best reference genes; 4re = all 4 
reference genes; bad2re = 2 least stable reference genes; EF-1a = EF-1o( as the reference gene; 18S = 18S rRNA as the reference gene. 
doi:1 0.1 371 /journal.pone.0097038.g008 



with minor modifications. RNA was precipitated in isopropanol 
overnight at 4°C. Tlie 8-aphid pooled sample was frozen with 
liquid nitrogen and ground with a clean glass rod for total RNA 
extraction in 400 |J,1 TRIzol reagent. The concentration and 
purity of RNA samples were measured with a NanoDrop-2000 
spectrophotometer (Thermo Fisher Scientific, Roskilde, Den- 
mark). The integrity was checked using electrophoresis in 1.5% 
agarose gels and 1 x Tris-acetate EDTA (TAB). The gel was then 
photographed in UV light (Bio-Rad, Universal Hood 11, USA). 

RT-PCR of Candidate Reference and Target Genes 

The total RNA isolated from plant leaves and the aphid samples 
of all developmental stages was reverse transcribed into first-strand 
cDNA only using the reverse primer of each gene. The 
components in each 25-|J,l RT reaction system included 7 |il 
RNase-free ddHjO, 1 |a.l 5 |a.M reverse primer and 2 |J,1 total 
RNA. The mixture was briefly centrifuged and denatured at 70°C 



for 10 mill, and then combined with 3.5 |J.l RNase-free ddH20, 
5 10.1 dNTPs (2.5 HM, TaKaRa, Dalian, China), 5 Hl 5x RT 
reaction buffer, 1 [0.1 M-MLV Reverse Transcriptase (200 U/|J.l, 
Promega, Madison, WI, USA), 0.5 (il Recombinant RNase 
Inhibitor (40 U/|J,1, TaKara). All reagents in the tube were 
centrifuged shordy and incubated at 42°C for 1 h. The RT 
reaction was terminated at 95°C for 5 min. The final cDNA 
samples were kept at — 20°C. 

The PGR reaction mixture for each fu-st-strand cDNA 
comprised 17.8 |0l ddH-^O, 0.25 |Xl forward primer (5 ^M), 
0.25 Hi reverse primer (5 [iM), 10 x LA Taq Buffer II (Mg^"^ plus), 
2 10,1 dNTPs, and 0.2 ^1 LA Taq (TaKaRa). The PGR program for 
all reactions was set as follows: pre-denaturation at 94°C for 
4 min; 32 cycles of denaturation at 94°G for 1 min, annealing at 
56-70°G for 1 min 30 s, and elongation at 72°G for 1 min 30 s; 
and 10 min at 72°G for the final elongation. RT-PGR products 
were purified and ligated with pMD 1 8T simple vector (TaKaRa), 
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Figure 9. Mean relative titre of YDVs (± SE, /7=3) in viruliferous winged adults of Rhopalosiphum pad/ after different virus-feeding 
durations. Virus titre is illustrated by relative expression values of the CP (A, C and E) and RTD (B, D and F) gene of YDVs. Expression values for the CP 
or RTD gene at each duration were normalized with the reference gene(s) selected by GeNorm and then compared with a one-way ANOVA (p) among 
these durations with each normalization condition. A and B: WYDV-GPV; C and D: BYDV-PAV; E and F: BYDV-GAV. 2re = 2 best reference genes; 3re = 3 
best reference genes; 4re = all 4 reference genes; bad2re = 2 least stable reference genes; EF-1o( = EF-la as the reference gene; 18S = 18S rRNA as the 
reference gene. 

doi:1 0.1 371 /journal.pone.0097038.g009 
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Figure 10. Mean relative titre of YDVs (± SE, n=3) in viruliferous wingless adults of Rhopa/osiphum padi after different virus-feeding 
durations. Virus titre is illustrated by relative expression values of the CP (A, C and E) and RTD (B, D and F) gene of YDVs. Expression values for the CP 
or RTD gene at each duration were normalized with the reference gene(s) selected by GeNorm and then compared with a one-way ANOVA (p) among 
these durations with each normalization condition. A and B: WYDV-GPV; C and D: BYDV-PAV; E and F: BYDV-GAV. 2re = 2 best reference genes; 3re = 3 
best reference genes; 4re = all 4 reference genes; bad2re = 2 least stable reference genes; EF-1cz = EF-lcx as the reference gene; 18S = 18S rRNA as the 
reference gene. 
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and finally transformed into competent DH5ot E. coli cells 
(TIANGEN BIOTECH, Beijing, China) according to a standard 
transformation method [61]. Three to five positive clones were 
sequenced by Sangon Biotech (Shanghai) Co., Ltd. 

RT-qPCR Data Collection 

Two-step RT-qPCR was used to obtain quantification cycle 
values (Cq, terminology from MIQE [62]) for all genes of aU 
samples. First-stranded cDNA of each sample was synthesized 
according to the protocol of FastQuant RT Kit (with gDNase, 
TIANGEN BIOTECH). Genomic DNA (gDNA) in a 20 |ll RT 
reaction system was first removed by adding 2 ^1 5x gDNA 
Buffer, 1 |i,g total RNA, and appropriate RNase-Free ddH20, and 
incubating at 42°C for 3 min. Then 2 IXIIO X Fast RT Buffer, 2 Hl 
FQjRT Primer Mix (random hexamers and oligo-dT primers), 
1 |ll RT Enzyme Mix, and 5 |ll RNase-Frcx- ddH^O were added 
and the mixture held at 42' C for 15 min, then at 95°C for 3 min. 
The fmal cDNA was kept on ice or at -20°C for the qPCR. 
Subsequently, SYBR Green I dye-based qPCR was performed for 
aU cDNA samples. Each 20 \iX qPCR reaction system contained 
6.4 III RNase-free ddHgO, 10 (il 2x SuperReal PreMix Plus 
(SuperReal PreMix Plus kit [SYBR Green I, TIANGEN 
BIOTECH]), 0.6 Hi forward primer (10 \lM), 0.6 Hl reverse 
primer (10 |lM), 0.4 (il 50 x ROX Reference Dye, 2 III dUuted 
cDNA (1:30 diluted cDNA sample template) or ddH20 (as the no- 
template control). qPCRs of all genes were run on a 7500 real time 
PGR system (Applied Biosystems) with the same program: pre- 
denaturation at 95°C for 15 min, and 3-step amplification 
procedure of 40 cycles of 95°C for 10 s, 59.2°C for 32 s and 
72°C for 32 s; followed by the melting curve stage of 95°C for 
15 s, 60°C for 1 min, 95°C for 30 s and 60°C for 15 s. The 
primers used for the qPCR reactions are listed in Table 2. In each 
single run, the raw Cq values were automatically calculated by the 
7500 real time PGR system detection software. Three ri^plicate 
qPCRs were done for each cDNA dilution sample. The 
amplification efficiency (E) of each candidate reference gene and 
the endogenous genes were calculated using the Cq values for a 5- 
fold dilution series of the cDNA (1:5, 1:25, 1:125, 1:625, 1:3,125, 
and 1:15,625). For virus genes, lvalues were determined using a 
2-fold dilution series of cDNA (1:2, 1:4, 1:8, 1:16, 1:32, and 1:64). 



can also determine the optimal number of reference genes by 
pairwise variation, F, between two sequential normalization 
factors containing an increasing number of reference genes, whose 
default cut-oflF value of 0.15 has been accepted by many studies of 
finding reference genes [15,44,63]. When values are higher than 
0.15, representing a large variation between the two sequential 
normalization factors, then another reference gene needs to be 
included to calculate a dependable normalization factor [25] . The 
pairwise variation starts from the first two and three (I'2/3) most 
stable candidate reference genes, terminating with the addition of 
the most unstable one. As the number of reference genes increases, 
the best number («) of reference genes for calculating the 
normalization factor wUl be decided, because the value of Vn/n 
+ 1 being first lower than 0.15, which means the variation in 
normalization factors from between n and n +\ reference genes is 
so small that an additional reference gene does not need to be 
included for calculating the normalization factor [25]. 

The normalization factors produced by GeNorm for the 
selected reference genes were used to calculate the relative virus 
titre and the expression value for the endogenous genes (as 
described in the GeNorm manual). To reach a persuasive 
conclusion, we used the best single reference gene (EE- la for aU 
experimental groups), the worst single reference gene (ACTl for 
WYDV-GPV-wingless group, and 18S rRNA for the remaining 5 
groups), a 4-reference gene combination (4re), 3-reference gene 
combination (3re), 2 best reference genes (2re) and 2 worst 
reference gene combination (bad2re) (according to Figure 2B) to 
normalize the expression levels for ago-la and dcrl in each group 
(using the serial formulas below ). The means among treatments 
were compared by a one-way ANOVA. 

ACqtarget = minCqtarget - sampleCqjarget ( 1 ) 



etarget = £'^^'"^«^' (2) 



Sequence Analysis 

The nucleotide sequences of all genes that were cloned from the . 
cDNA were analyzed using the BLAST tool on the NCBI website ^^^re = minCqre - sampleCq^e (3) 

(http:/ /blast.ncbi.nlm.nih.gov/) to confirm the genes. Vector NTI 
Advance 10.3 package (Invitrogen) served as the alignment tool. 



Algorithms and Statistical Analysis 

Expression stability of 4 candidate reference genes in YDV- 
viruliferous R. padi after the AAPs was evaluated by three 
algorithms: NormFinder 0.953 [37], BestKeeper version! [38], 
and GeNorm3.5 [25] respectively. Cq values were analyzed 
directiy with BestKeeper, and then the stability of the HKGs was 
ordered according to standard deviation (SD) values. In addition, 
Unear scale expression quantities (Q) were transformed from the 
Cq values by the delta-Cq method (ACq = min Cq — sample Cq; 
min Cq = lowest Cq value of the data set, sample Cq = Crj values 
of other samples) and corrected by the modified amplification 
efficiency [E = E +1) for the individual gene, according to the 
formula Qj^ E! ^'^'^^ then analyzed by NormFinder and GeNorm. 

To reach a reliable and unambiguous conclusion about qPCR 
results, we need two or more reference genes that are stably, i.e., 
consistentiy, expressed to normalize calculations [25,62]. In 
addition to ranking the expression stability of a gene, GeNorm 



6nre — V^2rer2re2^2re3^^^^^2r< 



(4) 



(5) 



RQ= 



Starget 



(6) 
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E'=E+l (7) 



All statistical analyses were done with the program SPSS version 
19 (IBM SPSS Statistics 19 Core System, USA). Plots were 
graphed using SigmaPlot version 12.5 (Systat, San Jose, CA, 
USA). 

Supporting Information 

Figure SI Total RNA extraction of oats and aphids. (A) 

Oat RNAs. H: healthy oat plants; BYDV-PAV: BYDV-PAV- 
infected oat plants; WYDV-GPV: WYDV-GPV-infected oat 
plants; BYDV-GAV-infected oat plants. (B) RNAs from mixed 
developmental stages oi Rhopalosiphum padi; (C) RNAs from BYDV- 
PAV- winged or wingless adult after various feeding durations (h); 
(D) RNAs from WYDV-GPV- winged or wingless adult after 
various feeding durations (h); (E) RNAs from BYDV-GAV- 
winged or wingless adult after various feeding durations (h). M: 
DL2000 (TaKaRa). 
(PDF) 

Figure S2 RT-PCR amplification of candidate reference 
genes and target genes. (A) RT-PCR amplification of the CP 
and RTD genes from total RNAs of BYDV-PAV-, WYDV-GPV- 
and BYDV-GAV-infected oat plants in Figure SI (A). H: healthy 
oat plants; W: ddH20 as no-template control; 1^: four YDV- 
infected oat plants; (B) RT-PCR amplification of 4 candidate 
reference genes from total RNA in Figure SI (B). Four lanes of 
each candidate gene indicated four RT-PCR products of 1/100 
diluted RNA sample under different annealing temperature; (C) 
RT-PCR amplification of endogenous genes from total RNA in 
Figure SI (B). Eight lanes of each gene indicated eight RT-PCR 
products of 1/10 (first four lanes) and 1/100 (last four lanes) 
diluted RNA sample under different annealing temperature. M: 
DL2000. 
(PDF) 
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